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EDITORS f PREFACE 


This translation has been prepared for the information of the 
atmospheric research community engaged in the application of general 
circulation models to the numerical simulation of climate. The model 
described herein has a number of features in common with other such 
models now in use. The parameterizations of short- and long-wave 
radiation transfer, the treatments of cloudiness and of convection, 
and the calculation of the fluxes in the surface boundary layer are of 
particular interest. Although no actual numerical solutions are given 
in this paper, the use of the so-called splitting method of Marchuk 
is envisaged. These and other features may make the paper of value 
to those engaged in similar modeling efforts. 

To the editors’ knowledge, this paper represents the most com¬ 
prehensive formulation of an atmospheric general circulation model 
undertaken in the Soviet Union as of its publication date (1970). 
Although the model is discussed here as a hemispheric one, the exten¬ 
sion to a global domain is straightforward. 




SUMMARY 


Numerical experiments, based on the physical laws of atmospheric 
motion and the processes of heat and moisture exchange, have attained 
importance for developing numerical methods of weather forecasting 

and for studying the genesis of climate and the general circulation. 

(1 2 ) 

Achievements in this field 9 unquestionably occupy an important 
place in the theory of forecasting and general atmospheric circulation* 
However, there is a basic inadequacy in the indicated studies, 
due to the fact that climatological values of moisture and cloudiness 
are used in icalculations of radiational heat fluxes. As a result, the 
dynamic characteristics and heat fluxes are not necessarily in balance. 

In this article is presented a scheme of the general atmospheric 
circulation that differs from the preceeding studies of the problem 
in physical content, since cloudiness and moisture are not given, but 
are determined in the process of integrating the appropriate equations. 
However, the ocean temperature is considered to be known from climate 
ological data. To reject assignment of this temperature would require 
that a system of equations describing the laws of oceanic circulation 
be included in the model. Consideration of such a complete problem is 
premature at present for a number of reasons. 

The general circulation in the northern hemisphere for the cold 
period of the year is considered in this study. As a result, it is 
possible to disregard evaporation from dry land, which greatly facil-t 
itates solution of the problem. 
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SYMBOLS 


The following basic symbols are used in the article, 

A - the flux of longwave radiation, downward 
a = the mean radius of the earth 
B = the flux of longwave radiation, upward 
C = the amount of cloudiness 
c^ = specific heat of air at constant pressure 
c v = specific heat of air at constant volume 
a = radiating capacity of underlying surface 
F - radiation balance 
f = relative humidity 
g = acceleration of gravity 
H = altitude of O surface 
H g = altitude of the earth’s surface 
L = latent heat of condensation 
2ft cos 0 « Coriolis parameter 
p = pressure 

P s = pressure at the earth’s surface 
Ttrp = turbulent vertical heat flux 
II q * turbulent vertical moisture flux 
Q = flux of total radiation 
q = specific humidity 
r « albedo of underlying surface 
= gas constant of the air 
S = flux of direct solar radiation 
T = air temperature 

T g = temperature of the earth’s surface 
t = time 

u = component of wind velocity along latitudinal circum¬ 
ference with positive direction to the east 
v = component of wind velocity along meridian with posi¬ 
tive direction to the South 
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A = longitude 

Y = coefficient of horizontal macroturbulent diffusion 
{also, scattering angle, Ed.] 
p = air density 


a 



coordinate along a vertical 


= component of frictional stress along latitudinal 
circumference 

Tq = component of frictional stress along meridian 
0 = co-latitude 

Q = angular velocity of earth's rotation 
dp 

w ^ = vertical-p velocity 
da 

^ = vertical-a velocity 

c 

K ~ ^ ” [also, scattering function, Ed.J 
v 


{Editors 1 note: There are also numerous other symbols defined 
locally in the text.J 
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I. BASIC EQUATIONS AND BOUNDARY CONDITIONS 
OF THE PROBLEM AND METHOD OF SOLUTION 


1. SYSTEM OF EQUATIONS 

Consider a system of thermohydrodynamic equations describing the 

transfer of momentum, heat, and moisture in the atmosphere. Let a 

global system of coordinates be used; the relationship of pressure, p, 

to its value at the earth’s surface, p , is given as an independent 

S (3) 

variable along the vertical (the so-called a-system of coordinates). 

The system of equations includes: 


1. Equation of motion 
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2. Equation of heat flux 
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3. Equation, of water vapor transfer 


■I 1 + — t + vV 2 q + £_ 

9t a sin 03^ a 30 do ^ 
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p 9a 5 
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(1.4) 


4. Equation of hydrostatic equilibrium, integrated along the 
vertical from a = 0 to 1 


.*?/ l d5; 


(1.5) 


5. Equation of continuity 
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Integrating the left part of this equation over the entire mass 
of the atmosphere and using the boundary condition 03 - 0 at O = 0 and 
a = 1, we obtain 


9p £ 

3t 


l 

/ 


Dd£ = 0, 


( 1 . 6 ) 


where for brevity 


D = 


a sm 


in 0 ^ 


9(p u) 9(p v sin 0) *| 
s + -- 


9X 


90 


(1.7) 


Finally, integrating the equation of continuity over the vertical 
from 0 to a and utilizing the known relationship 


dp s 

w = Pg S + 
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we obtain 


co 



DdC 



( 1 . 8 ) 


to 



MC + 0 f -H-j- + 2 

I a sin 0 dA ja 


3P £ 

90 


(1*9) 


Thus, considering that the components of frictional stress and 
the fluxes of heat and moisture appearing on the right-hand sides of 
the original equations are functions of the dependent variables and 
also of a series of given "external 11 parameters of the problem, we 
obtain a closed system of eight equations, (1.1) to (1.9) for deter¬ 
mining the eight unknown values u, v, T, q, p g , H, co, and to. We note 
that the sink of water vapor resulting from condensation and precipi^ 
tation, and the latent heat released in phase conversions of water are 
not calculated in the equations of water vapor transfer and heat flun* 
These processes are modeled in the scheme with the assistance of a 
special procedure, considered in paragraph 5. 


2. BOUNDARY CONDITIONS 

At the upper boundary of the atmosphere where 0 - 0 

(a) co = 0, (1*10) 


(b) A - 0, 

At the level of the underlying surface where a - 1 


( 1 . 11 ) 


(a) u = v = co = 0, 


( 1 . 12 ) 


For the ocean, the temperature of its surface is assigned 


T s - V X > ®>> 


( 1 . 13 ) 
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For dry land, the equation of heat balance for the underlying 
surface is used 


n T + 


lh + f = o. 

q 


Ci.14) 


It is assumed that flows of mass, heat and moisture reduce to 
zero at the equator, that is, where 0 = ^ 


v 0, 9Q 0, g0 0, ^ °* (1,15) 

All functions are assumed periodic with X. 

3. INTEGRATING THE EQUATIONS OF THE MODEL 

System Cl-1) through Cl*9) is solved by a method of separation. 

This method was first used by G. I. Marchuk for integrating equation® 
of atmospheric dynamics, and it was he who demonstrated its effective- 

ft) 

ness in a problem of short-term weather forecasting. In agreement 

with the principle of this method, a system of equations of general 
atmospheric circulation at each step in time is solved in several 
stages. The following is the sequence of the calculations. In the 
first stage, the changes in the meteorological elements caused by 
radiational heat flux, by vertical turbulent exchange and by inter¬ 
action of the atmosphere with the underlying surface are calculated. 

The values thus obtained are used as initial data for solving a system 
of advection equations describing the transfer of meteorological ele¬ 
ments along a trajectory of the flow, with horizontal large-scale dif¬ 
fusion included. Then, by solving a problem of adjustment, dynamic 
agreement of the resultant fields is produced. Finally, in the last 
stage, the cloudiness is determined, and the resultant fields of tern-* 
perature and moisture are corrected, with consideration of the released 
heat of condensation and of the processes of convection. 


4. SELECTION OF A GRID 

All variable values are determined at nodes of a latitudinal/ 
longitudinal grid with an increment AX along latitude and an increment 
A0 along longitude. 
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In applying a given model, AX = A8 = 5° is assumed, giving a 
total of 1297 nodes, including the pole, at a level in a hemisphere. 
Seven levels, corresponding to values of cr equal to 0.2, 0.3, 0.5, 0*7, 
0.85, 0.99 and 1.00 are used. The coordinates of each node are deter¬ 
mined by the index i, taken along a latitude to the east of the Green¬ 
wich meridian; by the index j, taken along a meridian in a southerly 
direction; and by the index k, taken along a vertical from the upper 
boundary of the atmosphere. Here, the step along the vertical is 
variable in size, depending on cf. Forecasting is accomplished with 
a time step At, with the index n determining the number of such steps. 
For the grid-point value of a function, the usual notation is introduced 


.12 


f i,j,k = 3^8, a k> nAt ^ 


where i = 0,1,...,1; j = 0,1,...,J; k = 1,2,...,K. 

The variable step along vertical Aa^. is equal to 


A \ = a k ~ a k-r 


a ,i6> 


In what follows, any of the indexes may be omitted when this does 
not cause incorrect interpretation of the equation. 




II. FLUXES OF RADIANT ENERGY 


It must be kept in mind, when developing a method for calculating 
radiational heat fluxes to be used in numerical models of the general 
atmospheric circulation and long-term weather forecasting, that hydro- 
dynamic methods give limited information on the state of the atmosphere. 
This is due both to the usually applied low vertical resolution in the 
atmosphere, and to the fact that certain parameters on which the radia¬ 
tion flux depends cannot be obtained from a hydrodynamic method. Of 
particular relevance are the characteristics of scattering processes 
in the atmosphere and the albedos of the underlying surface and of the 
clouds. Cloudiness has a special significance in calculating radia¬ 
tion fluxes. The means for calculating cloudiness now used in many 
studies demands knowledge of the degree of cloud cover. At the 

present time, this can be established only approximately, for example, 
by determining statistical relations between values of cloudiness 

obtained from observations and values of moisture or vertical velocity, 

(J 8) 

which are calculated in a numerical method. * 

Thus, it turns out that parameters known only approximately enter 
the calculation of radiation flux. In this connection, questions arise 
about the precision of calculation. The first: What sort of error Is 
introduced by inadequate knowledge of the magnitude of the radiation 
fluxes? The second: If the requirements for precision in calculating 
the fluxes of radiation are known, what precision is required when 
assigning other basic parameters? 

It follows that the precision necessary for calculating the fluxes 
of radiation depends on the sensitivity of the entire hydrodynamic model 
to variations in the flux of energy. 

An answer to the second question may be given only after numerical 
experiments with a hydrodynamic model including radiation heat fluxes. 

In view of certain individual experiments, the first of the ques¬ 
tions posed above can be studied independently of the hydrodynamic 
model. 
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1. FLUXES OF SHORTWAVE RADIATION 

The method for calculating radiation fluxes suggested by Manabe 

and Strickler is well known and is used. However, their plan has an 

essential inadequacy, since it disregards the transformation of the 

radiation field resulting from the absorption by aerosols and the 

effect of multiple scattering of shortwave radiation, 

C9) 

In an article a proposal is made to use the work of K, S. 

Shifrin, I. N. Minin, 0. A. Avaste and N. P. Piatovskaia^^ to 
construct a method for calculating the fluxes of shortwave radiation, 
taking into account the influence of multiple scattering. 

We shall consider that part of the spectrum from 0.29 to 0.72y. 

For a more complete consideration of factors acting in the atmosphere, 
we shall introduce absorption by ozone into the method for calculating 
the fluxes of total and direct radiation. 

Equation C2.1) is used for calculating fluxes to total radiation 
at the level z in the Interval 0.29 to 0.72y: 


0.72 

Q z = J J q CA)b*(A) 


-sec ilk,to ] 

e * Z dX. 


0.29 


( 2 . 1 ) 


Here J^dX is the intensity of solar radiation in the given part of 

the spectrum arriving at the upper boundary of the atmosphere; B*(A) 

is the function converting the flux to radiance at the level z 

(2 \ % ^ 

with repeated scattering v taken into account: 


B*(A) = 


2R (A)cos i 


4 + C3 - icpCL - r)x z (A) 


where 


R z = 1 + y cos i + (1 - y cos i I e 


CL » 

-sec i 

z 


Editors 1 note: Apparently reference [12] was Intended. 


( 2 . 2 ) 


(2.3) 
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The parameter characterizes the integration of the scattering 
function 


7r 


K, = 


f / K(y)sin y cos y dy, 

*4 


where y is the angle of scattering, K is the scattering function, and 
r is the albedo of the underlying surface. 

The optical density of the atmosphere is calculated according to 
the equation 


j 30 _ ^ \ -otz , x-n \ -Bz 

T Z T 0,p CX)e + T 0,a^ A 6 ’ 


(2.4) 


where Tq ^ is the Rayleigh spectral optical density, Tq ^ is the 
aerosol spectral optical density, a = 0.125 1/km, B = 0.898 1/km. 


CIO) 


In Eq. (2.1), k. is the logarithmic coefficient of ozone absorp- 

A C13) 

tion according to Vigroux, u) is the total content of ozone from 

z 

the upper boundary of the atmosphere to a given level z, according 
to Elterman (model 1964),^"^ and i is the zenith angle of the sun. 

The flux of direct solar radiation, S, through the horizontal 
surface, taking into account single scattering and the absorption by 
ozone in the spectral interval 0.29 to 0.72]i, is calculated according 
to the equation 


0.72 


/ -sec i[T (X)+k^o) ] 

Jq(X)cos ie dX, 


(2.5) 


0.29 


The fluxes of total and direct radiation at level z in the spectral 
interval 0.72 to 5.0]i are calculated by the equations 


5.0 


Q z - / J 0 O)B * O,T(5! H,0 )T(X C0, ) dX > 


0.72 
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- / J 0 


(X)cos ie 
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0.72 


t ( x h 2 o> t(x co 2 ) dA ‘ 


( 2 . 6 ) 


(2.7) 




The transmission function T(X), entering expressions (2.6) and (2.7), 
is calculated from the approximation equation 


T<X > - 1 - a(X)X + b(X) * 


( 2 . 8 ) 


Coefficients a and b for various parts of the spectrum for I^O and 
CC >2 are cited in the work of K. S. Shifrin and 0. A. Avaste. CIO) The 
calculation of the masses of the absorbing atmospheric constituents 
is made according to the equation 


\ = ^W k sec i(p + Z) n |' 


(2.9) 


where 


(p + it 




( 2 . 10 ) 



( 2 . 11 ) 


( 2 . 12 ) 


n is the coefficient expressing the influence of pressure variation 

on the value of the absorbing mass (n = 0.3 for ^0 and n = 0.4 for 

0° 

CC^); is the density of the absorbing material; is the content 
of absorbing material in a column of air from the upper boundary of 
the atmosphere to the level z, . 

Thus, with the help of expressions (2.9) and (2.10), the suspended 
absorbing masses of water vapor are calculated. 

The equations cited in Ref, 10 are used for calculating the ab¬ 
sorbing masses of 
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It is necessary to know the value of the flux of reflected radia¬ 
tion at various levels in the atmosphere. 

The intensity of the reflected flux is calculated from the equation 


BCe, Cp, i) = 



J o a)rCX)B*(X)T 1 (0, <p, i, X) dX. 


(2.13) 


For the flux of reflected radiation, we have 


R(i) 


2tt tt/2 

U 


B(0, cp, i)cos 0 sin 0 d0 dcp. 


(2.14) 


where 0 is the angle of reflection, cp is azimuth, and is a function 
describing the radiation absorption in the atmosphere. 

Disregarding the dependence on azimuth, in the spectral interval 
0,29 to 0.72y, the function is given in the form 


T i = 




sec i + 


(T^+k^oo^)sec 0] 


(2.15) 


In the interval 0.72 to 5.0y, the following expression was used 
for T^: 


TjCQ, i, X) = e 1 T [x^B, i, X)j T ^X^(9, i, X) 

The single-stage scattering of the ray reflected at the angle 0 
and the absorption by water vapor and carbon dioxide of the incident 
and reflected ray is calculated in Eq, (2.15). 

The absorbing masses and are calculated for the entire 

path traversed by the ray from the upper boundary of the atmosphere to 
level zk in accordance with the equation 


(2.16) 
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x' - x"(i) + x£(0). 


(2.17) 


0 ° 

where X 7 is the absorbing mass in the path of the ray from the upper 
boundary of the atmosphere to the reflecting surface Z, X^ is the 
absorbing mass in the path of the reflected ray from the reflecting 
surface Z to the level z^. 

The system of symbols used and the trajectory of the ray to the 

level z. are presented in Fig. 1* The function T(X f ) is calculated 

k k 

according to Eq. (2*8). The optical density t^ is computed by the 

equation 



(2.18) 


with the help of Eq. (2.4). 



Fig.1 —System of symbols 

In (2.14) we perform a numerical integration over 0. Here, follow- 
ing the study, } we separate the interval of integration into four 
parts: 0 to 30, 30 to 50, 50 to 70, 70 to 90°. Then, disregarding the 
dependence of opacity on azimuth, we obtain 


Editors- note: Apparently reference [9] was intended* 
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30° 50° 

R = 2ir B (20°, i) I cos 0 sin 0 d0 + B(40°, i) f cos 0 sin 0 d0 

' J 0 *4o° i 

(2.19) 

70° 90° . 

+ B(60°, i) j cos 0 sin 0 d0 + B(80°, i) j cos 0 sin 0 d01. 
*'50° */0° 


As a result, after calculating the integrals for the flux of re¬ 
flected radiation, we obtain the following expression: 

R = 2tt j B (20°, i)0.0670 + B(40°, i)0.1115 

+ B(60°, i)0.1500 + B(80°, i)0.1710 J . (2.20) 

We shall assume the albedo of the reflecting surface to be inde^ 
pendent of A. In this case, the intensity of the ray reflected at 
angle 0 is calculated by numerical integration over the spectrum. 

For 6 = 20° 

B(20°, i) = £ j J 0 C0.72)AXB*(0.72)T 1 (20 e , i, 0.72) 

+ J 0 (0.76)AXB*(0.76)T 1 (20°, i, 0.76) 

+...+ J 0 (4.5)AAB*(4.5)T 1 (20°, i, 4.5) j . (2.21) 

Expressions for B(40°, A); B(60°, A); B(80°, A) are written 
analogously. 

Aerological observations made by co-workers in the Department of 
Atmospheric Physics at Leningrad State University were used to verify 
the accuracy of the method proposed for calculating fluxes of shortwave 
radiation according to the suggestion of K. Ya. Kondratyev. Data 
from four flights were used for initial data. Observations were mada 
at 30 km above the earth’s surface. Three flights (H July 1964, 








11 July '64 


22 July *64 


23 Oct. '64 


21 Oct. ’65 



Measured 

Calculated 


— Fluxes of direct radiation through a horizontal surface 


11 July *64 22 July *64 


23 Oct. ‘64 


21 Oct. ’65 



- Measured 

-Calculated 

-R + 

-D, 


R, eal cm 2 min' 


Fig.4 — Fluxes of reflected radiation 
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Comparison of the curves shows satisfactory agreement between 
calculated and measured values. The greatest discrepancies are noted 
in the reflected flux values. This may be explained as follows. The 
flux of reflected radiation, consisting of the flux reflected by the 
earth’s surface and the flux scattered backward by a layer of haze ip 
the atmosphere between the earth’s surface and a given level, desig-j 
nated as D^, was measured in the course of observation. 

Equations (2.20) and (2.21) give only the first component of thp 
flux. Its second component is found if we calculate the haze by the 
equations suggested by I. N. Minin and K. S. Shifrin. However, those 
calculations are quite lengthy. We therefore used calculations made 
by N. P. Piatovskaia and K. S. Shifrin for different conditions 

of the underlying surface, and selected those which corresponded most 
closely to the albedo values of the underlying surface used in our 
calculations. Since opacity in haze depends significantly on albedo, 
approximation of this component may introduce errors in the value R. 

The shortwave balance, F^ = R - Q, in cloudless conditions cited 
in Table 1 was calculated from the values of total and reflected radia¬ 
tion. For comparison, the shortwave balance according to observational 
data is also presented. 

It is likewise interesting to compare the values of the radiation 
flux convergence in the layers. 

Profiles of the flux convergence, -AF^/Az, are presented in Fig. 

5, where AF.. was calculated as the difference in values of shortwave 
1 & 

balance at two levels. The value of -AF^/Az is identified with the 
middle of the layer. 

For quantitative evaluation of the discrepancies in the fluxes, 

the ratios S - /S and Q - /Q (Table 2) and the ratios 

calc. meas. x calc. meas. 

F. - /F- and AF- . /AF- (Table 3) were determined. 

1 calc. 1 meas. 1 calc. 1 meas. 

The negative value of flux is without physical meaning. It was 
obtained as a result of the nonsynchronous observations at different 
levels, connected with a change in the height of the sun. 


Editors’ note: The values plotted in Fig. 5 are not in complete 
agreement with the data of Table 1. 
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Table 1 

MEASURED AND CALCULATED VALUES OF SHORTWAVE RADIATION 



22 July 1964 

11 July 1967 

23 Oct. 1964 

21 Oct. 196i 


z, km 

F 

1 Cal 

F 1 Meas 

F 

1 Cal 

F 1 Meas 

F 

1 Cal 

F 1 Meas 

F 

1 Cal 

F 1 Meas 

0 

-0.3155 

-0.2830 

-0.3920 

-0.357 

-0.1850 

-0.1850 

-0.5347 

-0.5710 

2.200 

• ♦ • 

• • » 

• * * 

• * • 

#. • * 

• • • 

-0.5830 

-0.5910 

2.500 

» • • 

f • • 

-0.4732 

-0.4560 

• • • 

• • • 

1 « * 

« 1 • 


2.600 

-0.3828 

-0.4200 



» • i 

• * • 

• # • 

• • • 


4.000 

• • * 

• • » 

• • • 

• • • 

-0.2707 

-0.2610 

-0.6073 

-0.6090 

5.000 

-0.4770 

-0.4860 

-0.5437 

-0.5290 

... 

1 • • 

» « • 

• • • 



8.400 . ... -0.3293 -0.3520 -0.6278 -0.6480 

8.660 ... ... —0.6273—0.6020 ... ... ... • • • \ 



























Table 2 

COMPARISON OF CALCULATED AND MEASURED VALUES OF FLUXES 


21 Oct. 1965 


z, km 

S n 
cal 

Q cal 

s 

meas 

^meas 

0 

0.87 

0.942 

1.0 

0.97 

1.00 

1.6 

1.00 

1.01 

2.2 

1.02 

1.02 

3.0 

1.02 

1.02 

4.0 

1.04 

1.02 

5.0 

1.04 

1.02 

7.0 

1.07 

1.03 

8.4 

1.05 

1.02 

10.0 

1.05 

1.03 

14.0 

1.05 

1.05 

26.0 


... 

30.0 

| 

i 

i 

• • • 

... 


23 Oct. 1964 


z, km 

S cal 

Q cal 

s 

meas 

Q 

meas 

0 

• • • 

* ft * 

1.0 

♦ • • 

• • • 

1.6 

• • • 

ft • ft 

2.2 


0.99 

2.6 

1.08 

0.98 

3.0 

1.07 

0.97 

3.3 

1.04 

0.95 

4.0 

1.04 

0.97 

5.0 

1.00 

0.96 

7.0 

1.01 

0.98 

8.4 

1.01 

0.97 

10.0 

1.00 

0.97 

16.0 

0.99 

0.98 

20.0 

0.99 

1.00 


22 July 1964 

z, km 

S cal 

Q cal 

s 

%ieas 


meas 

0 

1.00 

1.05 

1.3 

1.14 

0.95 

2.0 

1.15 

0.98 

2.6 

1.16 

0.98 

3.6 

1.16 

0.98 

5.0 

1.12 

0.97 

7.0 

1.08 

0.97 

9.0 

1.06 

0.97 

10.0 

1.05 

0.97 

16.0 

1.00 

0.97 

20.0 

1.02 

1 

0.97 


11 July 1964 


z, km 

^cal 

Q cal 

s 

meas 

Q 

meas 

0 

1.24 

1.08 

1.3 

1.28 

1.06 

2.0 

1.19 

1.05 

2.5 

1.14 

1.03 

3.5 

1.11 

1.01 

4.0 

1.07 

0.98 

5.0 

1.10 

1.00 

6.0 

1.06 

0.98 

8.0 

1.05 

0.99 

8.6 

1.06 

1.01 

10.0 

1.05 

1.00 

11.0 

1.03 

0.98 

15.0 

1.01 

0.98 

27.0 

• * • 

1.01 

30.0 

• • * 

1.03 


Table 3 


COMPARISON OF CALCULATED AND MEASURED VALUES OF FLUXES 



23 

Oct. 1964 

z, km 

F cal 

AF .. 
cal 

F 

meas 

AF 

meas 

2.2 

1.0 

ft ft ft 

1.1277 

4.0 

1.0371 

ft ft ft 

0.6482 

8.4 

0.9355 

0.3969 

10.0 

0.8984 

ft ft ft 

• • • 

... 

• • ft 



11 July 1964 

z, km 

F cal 

AF , 
cal 

F 

AF 


meas 

meas 

0 

1.0980 

... 

0.818 

2.50 

1.0377 

... 

0.965 

5.00 

1.0278 

1.14 

8.60 

1.0420 

0.960 

10.00 

1.0384 

• * • 
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Table 2 shows that the discrepancy in values of fluxes is within 
the limits of 0 to 20 percent. At one point, the discrepancy amounts 
to 28 percent. The greatest differences are in the lower layers of the 
atmosphere up to 2 to 2.5 km. Apparently this circumstance may be 
partially explained by the disregard of the radiation absorption by 
aerosols in the calculation method. 

Although the problem of aerosol absorption has not been adequately 
investigated up to the present, indirect indications of its influence 
do exist. xhe problem has not received careful research. 

From Table 3, it follows that discrepancies in the values of flux 
can be quite significant. Analysis has shown that large errors are 
found in cases of small values of radiation heat fluxes. 

Thus, comparison shows that the present method for calculating 
the fluxes, S, Q, and R of shortwave radiation may be used for calcu¬ 
lation on observational data in actual individual cases and in a model 
of the general atmospheric circulation. 

The equations cited above were used to explain the influence of 
the underlying surface albedo, the scattering function, and the averag¬ 
ing of the zenith angle, and the effect of limited information about 
moisture on the value of the radiation flux in the spectral interval 
0.7 to 5y. 

For variations in the albedo of the earth’s surface, values of 
radiation flux with albedo values 0.16, 0.25, 0.6, and 0.8 are con¬ 
sidered in a numerical experiment. 

Fluxes of total and reflected radiation, effective flux F^ at a 
given level z, the flux convergence in the layer AF^, and the relation¬ 
ship of fluxes with different albedo values are presented in Table 4. 
From the table, it follows that the influence of the albedo of the 
underlying surface on the value of radiation flux convergence in the 
individual layers of the atmosphere is not great. Basically, the 
albedo influence is strongest in the lower atmospheric levels. Examples 
in the table relating to individual actual cases of water vapor distri¬ 
bution shows that a change of albedo from 0.16 to 0.8 leads to a dif¬ 
ference of flux convergence in the lower atmospheric level of a total 
of 25 percent. Thus, it follows that demands for precise assignment 
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Table 4 

COMPARISON OF FLUXES FOR DIFFERENT ALBEDOS 


<4-1 

o 











h B 






Albedo 




TJ * 

O M 

PQ C 

! 



0. 

,16 



0. 

,25 


M T3 

0) C 

P. P 

CL O 
& W 

Level 

Layer km 

Q 

R 

F i 

AF i 

Q 

R 

F i 

AF X 

10 , 

0.0 

0-2.5 

0,2433 

0.03918 

-0.2041 

-0.0795 

0.2445 

0.06112 

-0.1834 

-0.0818 


2.5 

2.5-5.0 

0.3164 

0.03276 

-0.2836 

-0.0547 

0.3166 

0.04873 

-0.2652 

-0.0553 


5.0 

5.0-8.6 

0.3700 

0.03164 

-0.3383 

-0.0605 

0.7202 

0.04657 

-0.3205 

-0.0606 


8.6 

8.6-10.0 

0.4302 

0.03136 

-0.3988 

-0.0179 

0.4303 

0.04612 

-0.3811 

-0.179 


10.0 


0.4481 

0.03132 

-0.4168 


0.4481 

0.04606 

-0.3989 


30 

0.0 

0-8 

0.2457 

0.03931 

-0.2064 

-0.1746 

0.2469 

0.06172 

-0.1852 

-0.1778 


8.0 

8-15 

0.4126 

0.03156 

-0.3810 

-0.0677 

0.4126 

0.04955 

-0.3630 

-0.0678 


15.0 

15-27 

0.4801 

0.03140 

-0.4487 

-0.1084 

0.4801 

0.04929 

-0.4308 

-0.1084 


27.0 

27-30 

0.5884 

0.03133 

-0.5571 

-0.0374 

0.5884 

0.04918 

-0.5392 

-0.0374 


30.0 


0.6258 

0.03132 

-0.5945 


0.6258 

0.04918 

-0.5766 



* 

3^ 
O M 

“5 

U T3 
<y c 
a 3 
p. o 

E> 03 


a 

i 

i 


Layer km Q 


Albedo 


0.6 


R 


AF- 


Q 



10 


0.0 


2.5 


5.0 

8.6 


0-2.5 


2.5-5.0 
5.0-8.6 


0.2490 0.1494 -0.09960 
0.3176 0.1258 -0.1958 
0.3706 0.1215 -0.2491 
0.4305 0.1204 -0.3101 


8.6-10.0 


30 


10.0 

0.0 

8.0 

15.0 

27.0 

30.0 


0-8 

8-15 

15-27 

27-30 


0.4483 0.1202 -0.3281 
0.2515 0.1509 -0.1006 
0.4129 0.1212 -0.2918 
0.4802 0.1205 -0.3597 
0.5884 0.1203 -0.4681 
0.6258 0.1202 -0.5056 


-0.0922 

-0.0578 

-0.0613 

-0.0179 


-0.1912 

-0.0679 

-0.1084 

-0.0375 


0.2518 0.2014 0.0504 
0.3182 0.1695 -0.1487 
0.3708 0.1637 -0.2071 
0.4307 0.1621 -0.2686 
0.4484 0.1621 -0.2863 
0.2542 0.2034 -0.0508 
0.4131 0.1633 -0.2498 
0.4803 0.1625 -0.3179 
0.5884 0.1621 -0.4263 
0.6259 0.1621 -0.4638 


-0.0988 

-0.0584 

-0.0622 

-0.0178 


1.029 
1.011 
1.002 
0.989 a 


1.160 

1.057 

1.008 

1.000 


1.243 

1.068 

1.017 

0.983 


-0.1990 

-0.0681 

-0.1084 

-0.0375 


1.017 
1.001 
1.000 
1.000 a 


1.093 

1.001 

1.000 

1.003 


1.138 

1.006 

1.000 

1.003 


Variations In values are connected with calculation errors 
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of albedo need not be severe for the calculation of radiation flux con¬ 
vergence, For example, selection of an albedo equal to 0.25 instead 

of 0.16 is permissible, as is a replacement of the albedo values 0.6 

(29) 

and 0.8 by their average value 0.7. 

Another important aspect of the methodology of calculating fluxes 
is the necessity of evaluating errors, caused by the scarcity of original 
data on moisture that are used for calculating the absorbing masses. Of 
the 32 known values that describe the profile of moisture, six were 
selected. The remaining 26 values, which were used for calculating the 
fluxes, were obtained by interpolation. 

Thus, radiation fluxes were calculated in two ways. In the first 
way—absorption by the water vapor mass was calculated according to the 
32 values of moisture taken from the observational data. In the second 
way—'the absorbing masses were found according to interpolated values. 

For interpolation, an exponential dependence of moisture on o was 

used, where O - p/p represents an approximate linear dependence of 

s 

pressure on height. 

The values obtained for the fluxes are presented in Table 5. Com¬ 
parison of these values shows that the use of interpolated values of 
moisture and pressure leads to a discrepancy in the values of the fluxes 
of direct and total radiation not exceeding several percent. Table 5 
contains results of calculations for two actual cases with a relatively 
smooth profile of moisture. Under actual conditions, singularities in 
the distribution of moisture usually occur with the appearance of clouds, 
in which case the absorbing water vapor masses must be calculated indi¬ 
vidually over and under the cloudy layers. 

Calculations with values of the parameter equal to 0.6, 1.2, 
and 2.0 were made to explain the influence of the form of the scattering 
function. Results show that the discrepancy in the value of total radia¬ 
tion amounts to 5 to 8 percent. 

The conclusions obtained in this section make it possible to use 
approximate values of the albedo of the underlying surface in calcula¬ 
tions of the flux of shortwave radiation, and to utilize a small number 
of levels in the vertical in calculating fluxes of radiation. 



Table 5 


OBSERVED AND INTERPOLATED VALUES OF RADIATION FLUXES 



Case 1 

Case 2 



S 


Q 


s 


Q 


1 

I 

r 

I 

i 

1 


1 


o 

4J O 

o 

+J o 

o 

u o 

o 

4J O 


fx a 

d a c 

cu a 

dad 

a a 

3 a. s 

a d 

3 p, (3 

1 — 1 

U 0 

o u o 

M O 

0^0 

0 

0 u o , 

u o 

o u o 

<1) 

j> 

*d a> «h 

*d a> *h 

*d a) *rl 



rd a) *h 

d a) -h 

,d a) -h 

<D 

£ d rt 

H f3 (U 

a 

•h a cd 

•h d co 

rl a td 

•h d rt 

•Had 

id 


£2 H H 

ts M 

[2 H H 

£5 M i—1 

£2 1 — 1 i — J 

[2 H rl 

12 H H 

0 

0.2886 

0.2864 

0.3451 

0,3426 

0.1664 

0.1653 

0.2124 

0.2110 

1.6 

0.3717 

0.3702 

0.3911 

0.3896 

0.2323 

0.2316 

0.2527 

0.2520 

2.6 

0.3882 

0.3959 

0.4089 

0.4067 

0.2750 

0.2705 

0.2847 

0.2800 

3.5 

0.4221 

0.4107 

0.4191 

0.4177 

0.2976 

0.2932 

0.3038 

0.2993 

5.0 

0.4264 

0.4258 

0.4307 

0.4301 

0.3238 

0.3223 

0.3279 

0.3264 

7.0 

0.4336 

i . ■ ' 

0.4343 

0.4365 

0.4372 

0.3511 

0.3491 

0.3539 

0.3519 

. . • | 

16.0 

0.4256 

0.4270 

0.4266 

0.4279 

0.4321 

0.4334 

0.4330 

0.4344 \ 

20.0 

0.4202 

0.4197 

0.4206 

0.4202 

0.4822 

0.4809 

0.4828 

0.4815 


2. FLUXES OF LONGWAVE RADIATION 

Fluxes of longwave radiation are calculated according to well 
known equations. 

For fluxes at level z, directed downward, under cloudless condi- 
tions, we use the expression 


E(M) 

A~ = E[TCm)J + / D(u - m)dE{T(u) ]. C2.22) 

E(jn) 

For fluxes directed upward 

E(m) 

B® = ElT(m)] - f D(m - u)dE[T(u)J 
0 

E(M) 

+ a - a) / 

0 


D(u + m)dE(T(u)], 


(2.23) 




- 22 - 


where 




(2.24) 


(2.25) 


and where D is the transmission function, E is the radiation of a per¬ 
fectly black body, m is the mass of the absorbing medium from the earth’s 
surface to level z, M is the amount of the absorbing material in the 
entire mass of the atmosphere, and is the density of the absorbing 
medium. 

The integral function, proposed by F. N. Shekhter, which considers 
the dependence of the absorption on the amount of water vapor, carbon 
dioxide, and ozone, is used for the transmission function. 

(18) 

A particular variant of Shekhter 1 s transmission function has 
the form 


D(H 2 0, C0 2 , 0 3 ) = D(H 2 0, C0 2 ) - A(0 3 ), C2.26). 

where D(H 2 0, C0 2 ) is the transmission function for water vapor and car¬ 
bon dioxide together, and A(0 3 ) is the absorption function for ozone. 

The expression for D(H 2 0, C0 2 ) is constructed so that the mass of 
water vapor is the only parameter, with the influence of C0 2 considered 
automatically on the basis of the functional connection between the 
masses of H 2 0 and C0 2 : 

D(H 2 0, C0 2 ) = D(W) = D*(W) - AD, (2.27) 


where 


D*- 0.471e-°- 695 ’ ¥ + 0.529e- 8 - 94 ^. 


(2.28) 
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Here W is the effective absorbing mass of water vapor, determined with 
the help of the expression 


z 

» - — f p v'p” dz. (2.29) 

The correct values AD are the following: 

W. 0.0001 0.0003 0.001 0.003 0.01 0.03 0.1 0.3 1 3 

100AD ... 7.4 8.2 8.5 8.9 9.2 7.4 5.2 2.2 0.2 0 


The calculation of the effective absorbing masses is performed 

(19) 

according to the equations obtained by R. L. Magrn from (2.29), 
using the equation of state and the hydrostatic equation adapted for 
calculation according to data on the water vapor and pressure at various 
levels: 


W. 

1 


i-1 


3=1 


(2.30) 


AW. = 
J 


0.01 


l. 

_1 


Vi 


) 


( Pj - 


? 3+l J 


(2.31) 


with i = 2,3,4,...; j - 1,2,3... Here i is the number of levels con¬ 
sidered in the calculation of absorbing masses. 

The absorption function for ozone is represented by the following 

data: 


"V ' ' 

0.01 

0.033 

0.06 

0.1 

100A(0 3 ) . 

0.4 

0.9 

1.6 

2.3 


3.3 4.1 4.4 
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Evaluation of the integrals in (2.22) and (2.23) according to 
Simpson's equation makes it possible to obtain approximate expressions 
for the fluxes.^ The effective fluxes of longwave radiation F 2 at 
each level are calculated in the form F 2 = B - A. 

3. CALCULATION OF FLUXES UNDER CLOUDY CONDITIONS 

”" _ (5 6) 

The method used for the consideration of cloudiness * consists 

of the following: the flux at level z is calculated as the sum of the 

fluxes relating to the cloudy and cloudless parts of the sky and is 

taken with corresponding weights. 

In three-layered cloudiness, fluxes from all cloudy layers enter 

the flux at level z, in correspondence with the scheme presented in 

Fig. 6 as an example. 



Fig.6 — Scheme of directed radiation flukes 
in cloudy conditions 


Generally, the flux from above, at level z, is calculated as the sum: 



+ 


(2.32) 
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Fluxes , directed upward from below are calculated according to 
the equation 



z 


(B S + R S ) 
z z 


x i=l ' 



+ ( B z 2 + R ! 2 ) e '2 + (\ 3 + E ! 3 ) e ' 3 - 


(2.33) 


The meaning of the quantities entering expressions (2.32) and 
(2.33) is presented below: 

^1 ^2 ^3 

A , A , A are fluxes of longwave radiation from the lower 
z z z 

boundary of clouds in the lower, middle, and upper strata respectively; 
k l k 2 k 3 

B , B , B are fluxes of longwave radiation reaching level z 
from the upper boundary of clouds in the lower, middle, and upper strata 
respectively; 

*1 ^2 "^3 

Q 7 , Q„ , Q are fluxes of shortwave radiation transmitted by 
z z z 

clouds in the lower, middle, and upper strata; 
k l k 2 k 3 

R , R , R are fluxes of shortwave radiation, reflected from 
z z z 

the upper boundary of clouds in the lower, middle, and upper strata 
respectively and reaching level z; 

A 

C_^ is the amount of cloudiness visible from below from given 


level z; 


2^ is the amount of cloudiness visible from above from given 


level z. 


The following relations are used for recalculating actual cloudiness 
visible from below or above, with i = 1 for clouds at the lower 
stratum, i = 2 for clouds of the middle stratum, and with i - 3 for 
clouds of the upper stratum. 
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C 


1 


C 


2 


A 



1 - Cc 2 + c 3 ) 



(2.34) 


(2.35) 


If one of the three cloud layers is absent, the corresponding 
value CL is assumed equal to zero and equations (2.32) and (2.35) re¬ 
duce to the appropriate particular case. 

We calculate the effective flux in cloudy conditions as 



4. CALCULATION OF RADIATIONAL HEAT FLUX 

To obtain an equation for calculating the fluxes with data at a 
limited number of layers, it is assumed that the effective flux may 
be presented as a quadratic function of the vertical coordinate in the 
form: 


F k = a + b (a - a k ) + c(a - a k ) 2 . 

We shall find expressions for the coefficients a, b, c in terms 
of the values of the fluxes at the levels k - 1, k 4* 1. If we intro¬ 
duce the designations 
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a k+l " a k A0 k+1 = p > 
a k ~ CT k-l = Aa k = q » 


we obtain 


b 



p 1 2 3 4 \-i - F k c q 2 - p2) 

pqCp + q) -* 


since 



and therefore 

3F I __ q _ p _ g - p 

3cT | cr=cr pCp + q) k+l q(p + q) k-1 ~ pq k‘ 
u k 

The net radiational heating (the radiational heat flux conver¬ 
gence) $ is expressed by the derivative 8F/9cr in the following manner: 
$ = 8/ c p^ s 9F/9cr. Fluxes are calculated from the following original 
data: 

1. Distribution of moisture, 

2. Distribution of pressure, 

3. Albedo of the underlying surface and of clouds, 

4. Absorbing properties of clouds. 
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III. INTEGRATION OF EQUATIONS OF THE MODEL 

1. CONSIDERATION OF THE VALUES OF FRICTION. HEAT, 

AND MOISTURE FLUXES 

Changes in the wind velocity, temperature, and moisture caused by 
vertical turbulent diffusion and the radiant flux of heat are described 
by the following equations 


3u _ _ j*_ 9T X 

3t p 3a * 
s 


(3.1) 


3v _ _ j*_0 

3t p 5ff ’ 
r s 


(3.2) 


9T = g 3_ 
3t c p 3a 

P s 


0I T + F), 


(3.3) 


a 3n 

= £ _ 1 

9t p 9a 
s 


(3.4) 


We shall assume that the height of the planetary boundary layer 
coincides with the level 0^2 > within the boundary layer, two 

calculating levels are located on and outside of this layer 

the intensity of friction and turbulent heat fluxes are equal to zero. 
Consequently, 


9u, 

k 

8v k n 

9q k 

9t 

3t °’ 

3t 


(3.5) 



9t 


P c 
s p 



- »(o k >- 


(3.6) 


Level C7 k coincides with the surface of the earth. 
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As a vertical coordinate, we shall take the geometric height z in 
place of a, and we shall write the basic equations for the boundary 
layer in the form 


3T 

9t 


3u 

3t 

1 OL X 
p 9 z 3 

(3.1') 

3u 

3t " " 

2 

u^ cos a. 

(3.10 

3v 

3t " 

2 . 

u* sin a. 

(3.20 

ku * T * ~ 

rrz [*00 - ta>], 

P P 

(3.30 

Ji_ . 

3t 

- ku*q*. 

(3.40 


Here a is the angle between the direction of the vector of the fric¬ 
tional force near the underlying surface the tangent to the circle 
of latitude. We also introduce the definition 


XCu, v, T, 



Xdz. 


Approximating x as a polynomial of the fourth degree in z, satis¬ 
fying all the boundary conditions (3.8), we obtain 

X = a 1 x(H) + a 2 x(h) + a 3 xCS), (3.12) 


where 


a 


1 


(6 - h)c 



a 2 (-5 - h) ( \ - c ) , 
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a 3 o-w||+c[i-(fE-|) 3 

6 - h 

C 20(H - h)* 

If CH - h) << 5, then with adequate precision we have 

a 1 = (6 ~ h)(c - 0.15), 
a 2 - 0$ - h)(0.25 - c), 
a 3 = 0.15 05 - h). 

Now, considering the boundary conditions (3.9) and the expressions 
(3.11) and (3.12), it is possible to write a system of equations (3.1 /7 ) 
to (3.4") individually for dry land and ocean in the following form: 

For dry land: 

8 2 

[a^u(H) + a 2 uCh)] = - u* cos a, (3.13) 

8 2 

[a^v(H) + a 2 v(h)] = - u* sin a, (3.14) 

ff [ ai T(H) + a 2 T(h)] = - - a 3 $(6). (3.15) 

P 

For the ocean: 

fjr [a^u(H) + a 2 u(h)] = - u* cos a, (3.16) 

[a^v(H) + a 2 v(h)] = - u* sin a, (3.17) 

ft £ a i T(H) + a 2 T(h)] = - ku*T* - [F(h) - F(6)] - a 3 $(6), (3.18) 

pP 
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[a-^H) + a 2 q(h)] = - ku^q*. (3.19) 

Since within the surface layer, the direction of the wind velocity 
changes little and coincides with the direction of the vector of fric^ 
tional force, it is possible to assume that 

cos a = u(H)/V(H), 

sin a = v(H)/V(H), (3.20) 

2 2 2 
where V = u + v . 

The system of equations (3.13) to (3.15) contains seven unknowns 

and the system (3.16) to (3.19) eleven unknowns. For completing these 

systems we shall use the known expressions for the profiles of wind 

( 21 ) 

velocity, temperature and moisture: 

u * 

u(h) = ■£-- ff(h/L) - fC£^/L]cos a, 

T(h) - T(z Q ) = T*[f(h/L) - fCz Q /L)J, 

qCh) - q(z 0 ) = q A [(h/L) - f(z 0 /L)]» (3.21) 

(here L is the length scale 

2 

u * 

L = ~ 2 --> 

(BT* + 0.61gq*) 

and 6 = g/T is the buoyancy parameter) and four empirical relationships, 
obtained in Ref. 22. 

• for ocean and dry land: 


- bu * T * 


10~ 2 ZV 2 (6) 

B 


yN(y), 


(3.22) 



~ 32 ^. 


• for dry land 


0 . 


u* = VOS)M(li) 0.089 - 0.00356 In 


V(6) 1 

Zz 0 J’ 


• for ocean 


u. = 0.030V(<5)M(y), 


,- 2 . 


-bu A q A = 10 V(.6) Iq (z Q ) - q(6)]N(y), 


Here 


N(y) = [0.274 + 0.0167y 1/3 ] 2 , 


M(y) 1 + O.SC^qq) °* 3 (i000^ + 0,08 ^1000^ ’ 


y = 6 


T(z Q ) - T(6) - 8.4 
ZV(6) ’ 


(3.23) 

(3.24) 

(3.25) 

(3.26) 

(3.27) 

(3.28) 


The empirical relationships (3.22) to (3.25) are not valid for 
lower latitudes. The question of the structure of the boundary layer 
near the equator still requires both theoretical and experimental 
studies. 


2. INTEGRATION OF THE ADVECTION EQUATIONS 

We shall write the system of advection equations in the following 

form: 


9u u 9u v 9ti - _9u 

9t a sin 0 9A a 90 W 9a 


9v ,_ u 9v v _9v - J9v 

9t a sin 0 9A a 99 W 9a 


1 9 2 u 1 9 2 u 1 

a 2 sin 2 0 9X 2 a 2 90 2 J ’ 

1 IJS l + 9 2 v ] 

a 2 sin 2 0 9A 2 a 2 90 2 J * 


(3.29) 


(3.30) 
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3T u 9T v 9T - 9T _ f 1 9^T 1 9 2 T 

9t a sin 0 9X a 96 W 3a 2 . 2 . „2 2„ 0 2 

La sm 0 3 a a 90 J 


l£L + _H_ 9a. + v9£-9a, = [_1_ 3^q 9^ 1 

9t a sin 6 3X a 36 3a 2 . 2 . ..2 2 _ Q 2 

L a sm 0 9 a a 90 J 


(3,31) 


(3.32) 


On the right-hand sides of these equations small terms of the form 

ptg L ff are omitted, since even for polar regions these components are 
a 

one or two orders less than the remaining terms* 

In turn, this system is divided into three sequentially solved prob¬ 
lems. First, values of u, v, T, and q are advected along the circle of 
latitude, then along the meridian, and finally along the vertical. For 
brevity of notation we shall introduce the so-called standard system of 
coordinates according to the relationships 


dx = a sin 0 dX, 
dy = ad0. 


(3.33) 


and the vector function 


u 

v 

i = T . 

q 

Then the system (3.26) to (3.29) assumes the form 

31 31 31 f 3 2 1 3 2 1" 

9t 3x 3y = V h x 2 + 9y2 _ * 


(3.34) 


(3.35) 


This vector equation is also solved by a method of splitting. It is 
divided into two equations in such a manner that derivatives are present 
only along one spatial variable. First we solve the equation 
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3 1 3 $ 3 2 1 



(3.36) 


during the time interval t <: t £ t + At, then, taking the function $ 
obtained as a result of the solution of (3.36) for the initial distri¬ 
bution, we solve the equation 


3 $ 3 $ 

’St’ V 9y 



(3.37) 


during the same interval in time t + At. Since the algorithm for solv- 

n 

ing these equations is the same, we shall limit ourselves to a consider¬ 
ation of (3.36). 

In schemes for numerically modeling the general circulation of the 
atmosphere and ocean given by Leithand Crowley, explicit finite- 

difference schemes of second-order accuracy along spatial coordinates 
for one scalar quasi-linear equation of the form (3.35) were proposed, 
and it was shown that these schemes are conditionally stable. However, 
with integration for a long period of time, they give significant errors 
in phase. Fromm^ and Crowley^ ^ developed more complex difference 
schemes, making it possible to lessen the errors in amplitude and phase 
of the numerical solution, as a result of increasing the order of approx¬ 
imation of the spatial derivatives. 

Unlike the referenced studies, we shall consider not one, but a 
system of equations (3.36) and find their common solution. In Ref. 27, 
a precise solution to the system (3.36) is found in integral form; by 
using it, it is possible to construct diverse finite-difference schemes. 

We introduce a new scalar function <J> and vector-function in 
agreement with the nonlinear relationships 

u — — 2v - ~ - n ¥-*•<(». (3.38) 

dX ~ ~ 

(27) 

It is not difficult to show, that d> and ¥ satisfy the equation of 

3Y 3 2, F 

heat conduction, for example ^ - V —? = 0. Considering the solution 

9t 3x 2 
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of this equation on an unbounded straight line during the interval of 

time t £ t £ t + At with initial conditions for (j) 
n n 


and for ¥ 



it/H iii in 

¥ = $ • <b , 

/v ^ T 7 


following from (3.38), we obtain for x = 0 


(3.39) 


(3.40) 


4> 


n+l 


1 




(3.41) 


^n+1 


1 

2/uvAt 


oo 







d?. 


(3.42) 


Hence, on the basis of (3.38), we obtain 



(3.43) 


Here, is the value of the function <J> at the origin. 

From (3.43) it is possible to obtain a solution at the nodes of 
the grid considered. We shall place the source at node i of the grid 
and assume there is a linear dependence of u on x in calculating the 
function (j> according to (3.39) 


~ c + bx. 


u 


(3.44) 
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We introduce the designation 


OO 

= f exp (- 4 ^)d?. 


(3,45) 


Evaluating this integral by parts, we obtain 


J k = 2vAt(k - l)J k _ 2 


°° 2 

-2v4t / e*p (3.46) 


1 n n 

On the basis of (3.38), substituting -g|— for - <j> n u n and using 
(3.44) and (3.45), we come to the following relationship: 


(1 + bAt)J k = - cAtJ" k ^ + 2vAt(k - l)J k _2< 


(3.47) 


Hence 


J, - - 3, 


J n = 


B 2 + y. 


J, = -3(6 + 3y), 


4 ? 2 

J A = 6 + 66 y + 3y • 


(3.48) 


Here 


- J k 

J k = 7^ 

k j q 


6 = 


cAt 


1 + bAt’ 


. _ 2vAt 
Y 1 + bAt‘ 


Now, if we introduce the function in the form of a portion of 
the Taylor series in the vicinity of x = 0; 


$ n - $ n = V* 


N 'aV^ k 
~ 1 x 
„ k / k!» 
k 3x /± 


k=l 
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and then substitute this expression in (3.43) and use (3.45), we obtain 



(3.49) 


From Eq. (3.49) a series of known approximations may be obtained 
as a particular case. We shall consider several such methods. 


1. c = u. , b = 0, N = 1. 
l 


In this case, in agreement with (3.49) 




(3.50) 


or 


a(^ - S^j^) with ^ 0 


n 


c 1 - < - 


- S’ 1 ?) with u n < 0 

'T+l ~L 1 


(3.51) 


a = 


u“ At 
~kT~' 


(3.52) 


What is obtained is an alternate scheme of first-order 
approximation. 

2. c = vu, b = 0, N = 2. 

Then 



(3.53) 


or, substituting derivatives of according to the equation 
of central differences and introducing the designation 


r 


vAt 

2 ’ 

(Ax; 


(3.54) 
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we obtain 



— ($ n 
2 l ~L+1 


$n J 

'^L-l 


+ (?- +r ) «5+i - 2 *i + Si-P- (3 - 55) 

(23) 

This equation was obtained differently in the work of Leith 
for the case v = 0. 

3. c = uj, b = 0, N = 4. 

In this case, approximating derivatives with fourth-order 
precision, we obtain 


$ n+1 _ $ n = _ £_ [8($ n - $ n ) - ($ n - $ n )] 

% ~L 12 iei ~L+l ~L-1-' ^~L+2 ~L-2- J 


- k + 2r > I 30 $S - 16( &1 + fi-l> + <£ + 2 + fi-2» 


12 C “ 3 + tor > I' 2 <£ + l - £-l> + fi+2 - Ji-2 ] 


+ \r (a 4 + 12a 2 r + 12r 3 ) [64? - 4(4?.. + 4? .) 
24 v ~l+ 1 'SL-l 


+ <*S« + S- 2 > 1 - 


(3.56) 


/on 

This scheme is presented in the work of Crowley for the 
case v = 0. 


,9u 


n 


Now we shall consider the case b = , c = vu. Assuming N = 2 

in Eq. (3,49), we obtain 


^n+1 -n 


- «1 ( $ n 
2 v Xi+l 


.n \ 

k-i* 


(£*o 


c$ n 

'■M+l 


2$ n + $ n ) 
~L ~L-1 ; 


( 3 . 57 ) 
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a' 


u" At 




Ax 


It is not difficult to show that this scheme satisfies the requirement 
of preserving the first moment of the equations of motion. 

We note that with the aid of (3.49), schemes can be obtained pos¬ 
sessing the property of preserving also the second moments for the equa¬ 
tions of motion. Thus considering 


3u . 3u 

9F +u K- 


- o. 


and assuming in (3.49) c = i (u. , - + u, + u. -) and b = 0, we obtain 

(33)* ^ 

the Arakawa scheme, satisfying the €%^ttions of preserving the 

2 

average values of u and u on a given with appropriate boundary 

conditions. 

We shall now consider the equation for vertical advection 


3$ _ 3$ 

9t = ~ w do’ (3.58) 

Taking 0) at the moment t = nAt, it is possible to obtain the equation, 
analogous to (3.55) 


,n+l ,n 

& ■ ik 


At(w£ At - Aa^) 

Aa k+l^ k+1 + Aa k) 


U n 
c Ak+1 



At(i \ &t + A< W ..a . .n 

Aa k(^ CT k+l + Aa k ) ^ 2k ~ 1 


Editors' note: See also A. Arakawa, "A Computational Design for 
Long-Term Numerical Integration of the Equations of Fluid Motion: Two- 
Dimensional Incompressible Flow, Part I," J. Comput. Phys 1(1): 119^ 
143, 1966. 
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where 


A0 k = a k " a k-l’ 

Aa k+1 = a k+l " V 

For a uniform grid (Acr^ = Aa^ + ^ = we obtain the formula 

(3.55) in which V = 0 and 

a? At 


3. SOLUTION OF A SYSTEM OF ADAPTATION EQUATIONS 

At this stage, the problem of the dynamic adjustment of meteor¬ 
ological fields is resolved. The following system of equations is 
considered 


9u ; = g 3H _ V 9p s 

9t v a sin 0 BX pa sin 0 9X 9 


^ ^ T R T 3p 

3v -7 g 3H B ^s 

at “ a 30 ' P a 36 ’ 

s 


H = 



(.3.59) 


(3.60) 


(3.61) 



(3.62) 


9T _ K - 1 TO) 

9t k p a 
s 


(3.63) 


Here and subsequently 


l = ('20 + 


u 

a sin 0 


j cos 


0 . 


(3.64) 
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The relationship (1.9) must also be added to the system for the deter¬ 
mination of the value of (a). 

We reconstruct the first two equations of the system. For brevity, 
writing 


f = 

u a sin 0 3A 


+ 


p a sin 9 9X * 
r s 


(3.65) 


we obtain 


f 

v 


_ 

a 90 


V 3 *b 

p a 90 ’ 

r s 



Iv 




(3.66) 


(3.67) 


If we introduce the complex velocity w = u + iv and the complex 
variable F = f^ + if^, the system (3.67) may be written in the form of 
one equation 


|| - iZw = - F. 


( 3 . 68 ) 


The solution of this equation has the form 


n+1 


n+1 n iZAt 

w = w e 


/ 


iZ(t 

Fe n+1 i%. 


(3.69) 


n 


Removing the average value of F from under the integral over the segment 
(t n , t ^) » which we designate by F*, we obtain 


n+1 


w 


n iZAt 

= w e 


iF* 

Z 


(1 - 


iZAtv 

e ) * 


(3.70) 
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Separating the real and imaginary parts of the expression (3*70) we 
come to the following equalities: 


n+l n n . 1 r * . V1 

u = u cos £ - v sm e - [f sm 8 - f (1 - cos £)], 

L U V 

v n+ l = u n sin e + v n cos e - \ [f (1 - cos e) + f sin e]. (3.71) 

l XX V 


where 


£ * ZAt. 


(3.72) 


« 5 * 

The values of p and T entering f and f are calculated at the time 
s v 

t + At/2 by the obvious scheme. Thus, we determine 


-L 

* n At f m,*. 

■ p . - ~ / d ds - 
J r\ 


(3.73) 


and 


n 


T* = T n + — K ~ 1 I- ifL 
2 ic p O 
r s 


C3.74) 


According to the determined values of T , we determine 


H* = H H-* 3 
s 


27 


dZ- 


(3.75) 


, * * * 

Then p , H and T are substituted in (3.65) and (3.66) and with the 
determined values of f and f and the velocity field at time t , the 
velocity components u n+ l and v n+ l are found from (3.71). Then 
P ^"sin 0 is determined according to the formula 


n+l n 


At 


s a sin 0 


/[ 


9 , * , 9 , * *■ . 

9A (u V + 96 Cv p s sin 0) 


dC. 


(3.76) 
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where 


* 


u 


* 


V 


1 , n+l , n N 

= j ( u + u ), 

1 f n+1 , n N 

= j C v + v ). 


From the values u n+ \ v n+ ^ and p^ + ^ we obtain a) n+ ^ from (1.9). In the 
same manner, the determination of the velocity field and surface pres¬ 
sure at the particular time is established. It remains to determine 
dynamically the adjusted temperature field at moment t^ + ^. This is 
calculated by the following equation 


T n +i = T n + (k - l ) . At T n (M n + ^ > C 3.7 7 ) 

ZKCTp 

S 

4. FINITE-DIFFERENCE APPROXIMATION OF DERIVATIVES 

In calculating derivatives along spatial coordinates, a constant 
step is introduced along longitude 


Ay = aA0 & 555 km 


(3.78) 


and a variable step along latitude 


Ax. = x.a sin 0.AX. 
3 J 3 


(3.79) 


Here x^ is an integral function which assumes the following values: 


x n = 12, x 9 = 6, x~ = 4, X/ = 3, x, =...= Xq = 2, 


x 9 X 18 


(3.80) 


Thus, the minimum value of the step Ax^ is equal to 390 km, and 
the maximum is 708 km. Derivatives along the horizontal are calculated 
with the aid of the central differences 



f l,J+l,k ~ f l,j-l,k 
2Ay 


f i+xj,j.,k ~ f i-xj ,j ,k 

2Ax. 

3 


(3.81) 



The value of D is approximated in the following manner: 




1 j (p s u) i +xj, j,k ~ ^s^i-x - j ^k 

2 ( Ax j 


, (p s v) i,i+l,k sin 9 1+1 ' (p s v) i.1-l.k sln 9 i-l 


Ay 


sin 9. 
3 


(3,82) 


Integrals appearing in the calculating formulas are computed by means 
of trapezoids. 

We now proceed to the determination of the desired functions at 
the pole, which is a singular point of the equations being considered. 
Assuming that the vertical component of vorticity and the velocity 
divergence are bounded at the pole^^ we obtain with 0=0 


u = uCcT, t)cos A - v(cf, t)sin A, 

v = u(cr, t)sin A + v(cr, t)cos A. (3.83) 

Here u and v are determined according to the values of the functions 
u and v on the circle of latitude closest to the pole, where these 
functions are calculated in the process of forecasting. For convenience 
of the calculations, the latitude $ = 0 is treated as the boundary of 
the area, containing 72 nodes. Components of the velocity vector are 
determined at these nodes by Eqs. (3.83), thanks to which the average 
value of the mass flow within the area at this boundary is equal to 
zero. Any scalar value is determined at the pole also by its value at 
the neighboring circle of latitude, and has one and the same value at 
all nodes of the boundary $ = 0. 
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iy. CALCULATING CONVECTION AND HEAT OF CONDENSATION. 
ESTABLISHING THE PARAMETERS OF CLOUDINESS 


Considered in this paragraph is a method for establishing the 
amount of cloudiness Cin tenths) and its boundaries, as well as a 
method for correcting the temperature, T, and the humidity, q, result¬ 
ing from the process of condensation and convection. 


1. FORECASTING CLOUDINESS 

Cloudiness in three strata is forecast. Here it is assumed that 
cloudiness in the upper stratum is in the layer 500 to 300 mb, cloudi¬ 
ness in the middle stratum is in the layer 850 to 500 mb, and cloudi¬ 
ness in the lower stratum is in the layer 990 to 850 mb. We will denote 
the tenths of cloudiness in the upper stratum C^, in the middle C^, and 

in the lower C, . The tenths of cloudiness in each stratum are deter- 
n 

A A A 

mined by values of the average relative humidity and f^, f 2 and £3 in 
the respective layer on the basis of the empirical graph of Smagorinsky, 
establishing the following relations: 


with 


f ± £ 0.60 


3.25 f, - 1.95 with 0.60 < f, <; 0.90 


with 0.90 < f- 


0.7 


with 


f 2 £ 0.35 


with 0.35 < f 2 ^0.85 

/A 

with 0.85 < f 0 


with 


f 3 £ 0.25 


C R = <'1.72 f 3 - 0.43 with 0.25 < f 3 <; 0.83 

with 0.83 < f 0 


(4.1) 


(31) 
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We shall tew? the values f ± =0.90, f 2 = 0.85, f 3 = 0.83 critical 
ones and designate them f^ cr » ^cr’ ^3cr* Average values of the rela¬ 
tive humidity in the layer are determined through the prognostic values 
of relative humidity according to the formulas: 


£, = 


f„ = 


f„ - 


(f, + f,), 

1/2 

2 \7 


h + h + 


J «5 + f 6>' 




(4.2) 


Here f^ are known functions of q^, T^: 


q k 

f = —-H 

k max 


v (T k ) 


max 


m HlCi A f 1 A 

where q^ (T^.) is determined by the Magnus formula 


max 3.8 ( T k " 273 

91, CT,,) zrz- ex P { 1*71 


k P s a k 


T. - 38 
k 


Processing to a determination of the boundary of cloudiness, we 
note that the vertical resolution assumed in the given model is in¬ 
sufficient to obtain actual values of cloud density with adequate 
approximation. Therefore, intermediate levels are introduced 


a k+l/2 2 ^k + a W-P 


k - 2,3,...,K-2, 


for which T, q and f are found by linear interpolation. 


% 

Editors' note: The intended expression is 


o o I T. - 273 

max ,„ s 3.8 k 

9,, CO = e xp 17.27 


’k Ps a k 


T k - 35.86 
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The following scheme for establishing the boundaries of the cloud¬ 
iness is adopted. If it turns out that there is cloudiness in any 
stratum as a result of calculations from Eqs. (4.1) and (4.2), then: 

a. In the event of unbroken cloudiness (C = 1), cloudiness is 
found on all successively situated layers (primary and inter¬ 
mediary) corresponding to the given stratum, for which 

f > f cr beginning at the level with maximum f; 

b. In the event of broken cloudiness (C < 1), cloudiness is 
ascribed only to that level where f is a maximum. 

From each level defined in such a manner, cloudiness spreads up¬ 
ward and downward to 1/4 the thickness of the primary layer, that is, 
the layer included between the two neighboring basic O levels. Possible 
variations in the disposition of cloudiness are presented in Fig. 7. 
Primary levels between which cloudiness is found are represented by 
heavy lines. Thin lines indicate boundaries of cloudiness. An area 
occupied by cloudiness is cross-hatched, while cloudiness of various 
strata within the limits of one primary layer is indicated by cross- 
hatching with a different slant. 

In this scheme, all levels were considered downward from above 
(from 0 2 to 0^^). If cloudiness of a particular stratum is noted at 
a given level, that level is no longer considered in the distribution 
of cloudiness in the lower-lying stratum, which excludes the possi¬ 
bility of the presence of cloudiness at two strata at the same time. 

The scheme constructed for prognosis of cloudiness is completed 
by the following conditions: 

T max ^0.12 . 

a. in case q. < ——, cloudiness at level k is absent: 

k p 3 a k 

b. At level K - 1 cloudiness is always absent. 

2. CONSIDERATION OF CONDENSATION AND CONVECTION 

In the present scheme, processes of condensation and convection 
are considered individually in two stages. In the first stage, a 
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Fig .7— Possible distribution of cloudiness along tbe 
vertical in a primary layer 


correction made of the temperature and humidity fields as caused only by 
condensation in cloud layers; in the second stage, a correction is made 
as caused by convection, where it is assumed that the specific humidity 
in the cloud during the entire process of convection remains saturated. 

Consideration of condensation . Condensation in layers of cloudi¬ 
ness lying in primary a levels is considered. If the value of the 
relative humidity at level a^, containing cloudiness of the \x stratum, 

exceeds f , the heat of condensation is released, as a result of 
per* 5 

which the temperature T v and the humidity q, of the air change. Their 
new values and q^. are calculated from the formulas: 



max 

q k 


C V + 6 k CT 


* 

k 



(4.3) 



XG V C <k ‘ q k } ’ 


(4.4) 


where 



dq rx> 

dT 


X 


c 

P 


Equation (4.3) is obtained from the condensation condition 


q k 


f max (t ) 
ycr M k v k 


(4.5) 
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maX 7? 

If (T^) are separated into a series by stages T^, - T^. and are 
limited by two primary members of the separation, which in the given 
case insures sufficient precision, then we will obtain Eq. (4.3)* 

Equation (4.4) is obtained in accordance with the suggestion of 
Smagorinsky that the process of condensation occurs so rapidly that 
only a local change in temperature takes place. Then the energy equa¬ 
tion may be written in the form 

C p If At = “ LC y(9* - q)- C4.6) 

* 

Substituting the increment of temperature with - T^, we obtain 
Eq. (4.4). 

Unlike existing schemes for calculating the latent heat of con¬ 
densation, here it is foreseen that condensation may occur even with 
broken cloudiness. 

Calculation of convection . For calculating convection, a correc- 

tion to the temperature and humidity are made if the vertical gradient 

of temperature in saturated air exceeds the moist adiabatic , or in 

ad 

unsaturated air the dry adiabatic y . In each primary layer, the 
values y 11 ^ and y ac ^, calculated in the a system, are ass tuned constant 

along the vertical within the limits of this layer. Temperature values 

• xn a ad 

entering the expressions for y and y , are selected as averages, 

according to the corresponding primary layer. 

In order that the influence of any convective process may be ex¬ 
pressed every time by a change in temperature and moisture in primary 
levels (since only the values T and q may be considered in other parts 
of the scheme) it is assumed that the temperature changes in proportion 
to y ma in the cloud and in proportion to y ac * outside the cloud in each 
primary layer. In this case, if convection occurs, it completely en¬ 
velops the primary layer or several such layers. We designate the 
amount of cloudiness occupying the upper, middle or lower part of the 
parimary layer by C^, and C^, respectively. Then, the condition 

for the absence of convection may be written as: 

# 
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k+1 


- T, 


A a. 


ad,. 

Yv (1 " 


C) + 


ma^ 
Y k C, 


(4.7) 


where 

C = i-CC^ + 2C 2 + c 3 ). 


Convection occurs in the event of non-fulfillment of the condition 

(4.7), as a result of which values of the temperature, T, and moisture, 

q, in the primary layers change. If convection envelops m primary 

& 

layers, then new values of temperature T at corresponding primary 
levels may be defined from the system of equations: 

c±) + Y^ a C i> , i = k,..., k + m-1, (4.8) 


* * 

T - T 

-i±I_i = Y ad f _ 

Aa. ^ 


E 

-i=V 


( T J 4- I* +1 )A ai 


=E 


(T ^ ) Acr^. 


(4.9) 


Equation (4.9) is obtained on the assumption that convection occurs 


without condensation. In this case 


(32) 


the condition is fulfilled: 



(T - T )dCT = 0. 


(4.10) 


Evaluating this integral approximately by the trapezoidal equation, we 
obtain Eq. (4.9). 
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